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Abstract 

In this paper we present a treatment of hypersingular integral equa- 
tions, which have relevant applications in many problems of wave dynam- 
ics, elasticity and fluid mechanics with mixed boundary conditions. The 
main goal of the present work is the development of an efficient direct nu- 
merical collocation method. The paper is completed with two examples 
taken from crack theory and acoustics: the study of a single crack in a 
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linear isotropic elastic medium, and diffraction of a plane acoustic wave by 
a thin rigid screen. 



1 Introduction 

Many problems of fluid mechanics, elasticity, and wave dynamics (acoustics) with 
mixed boundary conditions can be reduced to hypersingular equations of the 
following type 
1 



f<P(t) 

-l 



- j + K o OM) 



(x-t) 



dt = f(x), \x\<l, (1.1) 



where K (x, i) is a regular ( in some sense) part of the kernel . Direct numerical 
treatment of eq.(l.l) is not easy. Many authors applied the approach, which may 
be currently considered as classical. If one defines a solution bounded on the 



both endpoints x — ±1 ( that will be shown below to vanish as y\x± 1| at 
x — > ±1 ), then one may search it in the form of the series 

oo 

ip (x) = Vl - x 2 ^2 ¥j Uj (x) , \x\ < 1 , (1.2) 

3=0 

where Uj (x) are the Chebyshev's polynomials of the second kind [1]. Then, using 
the relation 



i 



/ Vl t l u j (t)dt = -4- f Vl 12 U, (t) dt = 

J (x-t) 2 1 dx J x-t 

= -7r-^[T j+ i(z)] = -7rO- + l) U 3 (x) 



( Tj(x) is the Chebyshev's polynomial of the first kind), applying the scalar 
product of eq.(l.l) with the functions \J\ — x 2 Ui(x), one can reduce initial 
integral equation (1.1) to the infinite linear algebraic system of the second kind: 

7T 2 °° 

+ + S k a Vi = /i > i = 0,l,..., (1.3) 

Z j=0 



kij = / / Vl-x 2 Vl-t 2 Ut (x) Uj (t) K (x, t) dxdt . (1.4) 



-i 



We have used here the well known orthogonality relation [1] 
1 

J Vl - x 2 Ui (x) U j (x) = ^5 ij , i,j = 0, 1,-, 



-i 



where 5ij is the Dirac delta. 



The study of qualitative properties of the infinite system (1.3) is not the aim 
of the present work, and we only note that, when the regular kernel K (x,t) is 
of a convolution type [K (x,t) = K (x — t)], this approach is often applied in a 
different way. 

We rewrite eq. (1.1) in equivalent form 

1 oo 

J K (x — t) ip (t) dt = f (x) , \x\ < 1 , K(x) = — J L(s) e~ isx ds, (1.5) 

— 1 — oo 

where L (s) is the Fourier transform of K (x). Let L (s) be integrable for finite s, 
and for s — > oo it has the following asymptotics 

L(s)= A\s\+L (s) , L (s)=o(^), (5 > 0) , (1.6) 

then 

2A 1 r°° 

K(x) = + K (x) , K (x) = — L (s) e~ lsx ds , 

X 271 J-oo 

where K Q (x) is at least continuous, and the generalized value of the integral [2] 
has been used, that is, 



oo 

s\ e~ lsx ds = 2 lim / e~ es s cos (sx) ds — - 

oo 



oo „ 

f - 2 
e es s cos (sx) ds = -. (1.7) 

J X 2 



Equation (1.5) together with (1.6) is evidently equivalent to eq.(l.l). If one 
substitutes the series expansion (1.2) into eqs.(1.5) and then applies a scalar 



product with functions \J\ — x 2 Ui(x), one finally arrives at the infinite linear 
algebraic system [1] 

f; k m = fi, (t = 0, 1, ...), l l3 = || ^(l-x 2 )(l-t 2 )U t (x)U 3 (t)dxdt 

(1.8) 



3=0 



oo 

X 

— oo 



/ L(s)e-^ds ~ / J - +1(S)J ^ 2 (S)L(S) ^ S . 



equivalent to system (1.3). Note that the last integral is finite for all i,j = 0, 1, 
since J n (s) ~ O (s n ) , s — > 0; J n (s) ~ 0{\/y/s), s — > oo, which, on taking 
into account the asymptotic property (1.6), is evident. 

The integral equation (1.1) with a hypersingular kernel can thus be reduced to 
an infinite system of linear algebraic equations. This implies numerical treatment 
of double integrals in (1.4) or single integrals over infinite interval of strongly oscil- 
lating functions, like in (1.8). When solving integral equations with more regular 
kernels, a powerful direct collocation technique, in the framwork of the Boundary 
Element Method [3], can successfully be applied. For the Cauchy-type integral 
equations, a direct collocation method has been developed by Belotserkovsky and 
Lifanov [4, 5]. Hence, the main goal of the present paper is to develop an efficient 
direct numerical collocation method. 



2 Some properties of hypersingular integrals 



First, one should clarify in which sense hypersingular integrals as given by eq.(l.l) 
may be treated, since they do not exist either as improper integrals of the first 
kind or as Cauchy-type singular integrals. At least three different definitions of 
hypersingular integrals are known (see [6]): 

1. The integral is a derivative of the Cauchy principal value 

<p(t) ^ d }<p(t) 



J (x — t) dx J x — 



-dt. 



(x — tY dx J x — t 

a y ' a 

2. The integral is treated as a Hadamar principal value [4,5] 



(2.1) 



I 



<p(t) 

(x - ty 



-dt = lim 

£^ + 



x+e . 



Lp (t) dt 2<p (x) 
(x - tf ~ 



(2.2) 



3. The integral is a residue value, in the sense of generalized functions [2], 
that is is an analytical continuation of the integral 

b 

\x - t\ a tp (t) dt (2.3) 



/ 
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from domain, where it exists in the classical sense, to the value a = —2 . 
When ip (x) = 1, the three different approaches give the same result: 



j dt d j dt d j / b — x\ a — b 

J (x — i) 2 dx J x — t dx \x — a J (x — a) (b — 



dt 



(x - ty 



lim 

e^+0 



f f f f 

+ - + - + 



x — a e e x — b 

b 



X) 

a — b 



(x — a) (b — x) 



J \x-t\ a dt = J (x-t) a dt + J (t-x) a dt = 



x — a) (b — x) 



a + l 



+ , (2-4) 

a+l a+l ' v ; 



when Re(a) > —1. The analytical continuation of (2.4) from the half-plane 
Re (a) > — 1 to the value a = — 2 gives 

b 

dt 



I 



lim 



j (x-ty —2 



\a+l 



(.r - (i) (b — x) 



a+l 



a + l 



a + l 



a — b 
(x — a) {b — x) 



All three definitions are equivalent to each other, if the density ip (x) is ana- 
lytical (i.e. differentiable) on the open interval (a, b). Really, let p(x) be analytic 
and Re (a) > 1, then 



d_ j up (t) dt _ _ ^ d_ /V j \ p (t) dt 
dx J x-t s^+odx\J J j x-t --+" 

a \ a x+e J 



lim 



p{t) dt 



' x~e b 



2p (x) 



lim 



' x—e b 



— lim 

dx 



, J \x-t\ a p (t) dt= lim o I J + J J | a; - t\ a p (t) dt 

a \ a x+e I 

b 

(x - t) a <p (t) dt + J (t-x) a p{t)dt 

x+e 

' ~ t) —p(t)dt+ / [t ~ X J p(t)dt 

1 J a + l 



a 



x+e 



so, by applying analytical continuation to the last relation, one can see (with 
the use of a standard (e,5) formalism) that the right-hand side results in (2.1). 
Therefore, equivalence of 1 to 2 and 3 is evident. Let us prove that if <p (x) e 
C 2 (a, b), then a finite value of the limit in expression (2.2) exists, and so for 
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x € (a, b) the integral J a b ip(t)/(x — t) 2 dt is finite in any sense. Really, expression 
in the square brackets in eq.(2.2) is 

which has a finite limit at e — > +0. 



3 Integral equation with characteristic hyper- 
singular kernel 

Consider hypersingular equation with the characteristic kernel 

jpt)ti =fixh xe (a,b), f(x)EC 2 (a,b) . (3.1) 

Let us prove that a bounded solution of Eq.(3.1) is unique and is given as follows 



x-a)(b-x) b r f(t)dt 



y {x -a)\v- x) r 
9(*) = -2 J 



J y/(t-a) (b-t) (x-t) ' 
Really, according to the definition (2.1), equation (3.1) is equivalent to 

a a 

where C is an arbitrary constant. Now an inversion formula for the Cauchy 
characteristic integral operator [7] determines the bounded solution as 



J( x -a)(b-x) } f(t) 

g (x) = * / =- -dt, x e (a, b) , 3.2 

^ •> . /(* - a) (b - t) (x - t) 



and the constant C as 



f 2 i ^ - a) (ft - i) 

Thus, as indicated in the Introduction, any bounded solution of Eq.(3.1) vanishes 
at x — > a, b. 
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To construct a direct collocation technique to solve eq.(3.1) for arbitrary right- 
hand side, we divide the interval (a, b) into n small equal subintervals of the length 
h — (b — a) /n, by the nodes a = to, t±, ti, ■ ■■,t n -i,t n — b , tj = a + jh, j = 
0, l...,n. The central points of each sub-interval (ti-i,ti) are denoted by Xi, thus 
Xi — a + (i — 1/2) h , i = l,...,n. 

If we try to arrange an approximation of integral in (3.1) by using finite sum, 
like in regular cases, then for x = X{ we have 

rb g(t)dt A > */ dt 



a ( Xi - ty 

1 



j=0 ti-1 



X i tj Xi tj—i 



(Xi - 1) 

n 



V2 dt 

-h/2 t 2 



(3.3) 



Xi - 



where a value of the above hypersingular integral of the function 1/t 2 has been 
used. So, we try to approximate eq.(3.1) by the linear algebraic system 



3=1 



K tj %i — 1 



l,...,n. 



(3.4) 



Further considerations are similar to those in [4,5]. Let us prove that, by 
assuming x = x\ e (a, b) is fixed, the difference between solution g{xi) of the 
system (3.4) and analytical solution (3.2) tends to zero, when n — > 00. Obviously, 
the system (3.4) can be rewritten as 



J29(t 3 : 



Xj t{ Xj ij— 1 / 



and its principal determnant is 
1 1 



X\ — t\ X\ — to 

1 1 



X\ t n X\ t n —\ 



r (xi) , 



l,...,n 



X n t 



n-1, 



(3.5) 



(3.6) 



Let us rewrite the i — th row (i = 2, ...,n) of A as a sum of all other rows, 
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from k — 1 to k — i: 
t\ — t 



ti — t 



x 



(Xi - ti) (Xi - t ) 

tn ~ to 
(Xl ~ t n ) (Xi ~ t ) 

1 1 



(x n - t ± ) (x n - t ) 
tn ~ to 

\X n tn) ^Xn 



n (ti - to) 

3=1 

n 

n (xi - 1 ) 

i=i 



X 



x x - ti 

1 



X n ti 



(3.7) 



n(tj-to) uu(t q -t p )uu(x P -x q 

j q<p q<p 



nn (x q - t p ) 

q, p 



II (xi - t ) 

Xi t n 

where the lower limit in all products is 1 and the upper is n. A known value of 
the last determinant has been taken from [4,5]. 

According to the Cramer's rule, one needs to calculate the determinant A; 
where a column with the right-hand side (3.5) is substituted instead of the / — th 
column of A (3.6): 



1 



1 



Xi — ti Xi — to 



Xi ~t 2 Xi - ti 

1 1 



\Xi t n Xi t n _i 

h-tp 

(xi - ti) (xi - t ) 

t2 — to 

(xi - t 2 ) (xi - t ) 

t n — tp 

(xi - t n ) (Xi - t ) 



... f'(Xi) ... 
... f'(x 2 ) ... 

... f'(Xn) 

f'(xi) 

f'(Xi)+f'(x 2 ) 
n 

E /' M 

k=l 



1 



l l 



^x n t n X n t n _i 

ti-to 



(x n 


- h) (x n 


-to) 


t 2 — to 




\X n 


- h) (X n 


-to) 




t n — to 




(x n 


tn) (x n 


-to) 



where we have used the same summation of rows as in the case of A. 

The last determinant may be calculated arranging expansion by elements of 
the / — th column as follows (see also [4,5]) 



A, 



II(f;-fo) 
J_ 

n (xi - 1 ) 



E 



E f(x k ) 



k=l 



-i t m to 



("I) 



n n (t q -t p ) n n (x P -x q ) 

m+l 1<P< q,P^m q<p; q,p^l 

n n (x q -t p ) 



so 

m 

A e /'(^fe) n (xi - t p ) n (x q - t m ) 

9{xi) = ^ = {xi ~ to) £ (tm - to) S - J) n (* m - 1) n S - ^) ' (3 ' 8) 

The last expression at /i — > ~ n — > oo, under the condition xi G (a, 6) is 
fixed, can be estimated as follows: 

n l—l n 

n (xi - t p ) n {xi - t p ) n (x t - t p ) 

P =\ _ , = {xi -f,)f? x ^ = 

* {Xq ~ Xl) n ^ (3 . 9) 

if x belongs to the open interval (a, b) (i.e. I ^ n, 1^1). Here we have used 
the asymptotic estimate [4,5] 



m=l 



Further, by analogy 

; r, f] = (-ir(t m -x m )u i+i) n i-i~(-ir^M. 

I! (t m -t p ) q=l V 2g/ g A = A i V 2g/ n^b - t m 

Other terms in eq. (3.8) can be simplified as follows (h — > 0): 

TO *"» 

(xj - t ) -»• (arj - a), (i m - t ) -> (t m - a), /i ^ /'(i fc ) -> / /'(i)cft = /(t m )» 



fc=i 



hence expression (3.8) with /i — > tends to 



fc 771 7^ 



9 (zj) ^(xi-a) (b-xi) 



71 m=l y/ (t m -a) (6 - *m) 0^-0 

- a) (b - Xl ) ) f( t )dt 



n 2 



J fi - a )(b- t) ( Xl - t) ' 

An example on application of the proposed numerical method is shown in 
Fig.l, where (a, b) = (—1, 1). 
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4 Full hypersingular equation 



Consider the full equation 

b 

i 

-2 + K (x, t) 



I 



(x-ty 



where 



K (x,t) 



g (t) dt = f (x) , x G (a, b) , 

dK x (x, t) 
dx 



(4.1a) 



(4.16) 



Its bounded solution can be constructed by applying inversion of the characteristic 
part, that reduces eq.(4.1a) to a second-kind Fredholm integral equation 





g(x)+ J Ni (x, t) g (t) dt = fr(x), x G (a, 6) , 



(4.2) 



where 



Wi(x,t) 



h (x) = 



\J (x — a) (b — x) 



7T Z 



K x (r, t) dr 



a V( r - a )( 6 - r ) ( X ~ T ) 



<J(x-a){b- x) 



f f(r)dr 
I y/(r - a) (6 - r) (x - r) ' 



Let us prove that, if f (x) G Ci(a,6); K\ (x,t) G C\ [(a, b) x (a, 6)], then for 

any x G (a, 6) the difference between solution g (x) of the linear algebraic system 
1 1 



E 



%i tj—1 



+ hK (xi,tj) 



g(tj) = f'{xi), i = l,...,n 



and the bounded solution of eq.(4.1) tends to zero when h — > (i.e. n — > oo). 
Really, if one transfers the terms, related to the regular kernel, to the right-hand 
side and solves so written linear algebraic system with the characteristic matrix 
1/ (xi — tj) — 1/ (xi — tj-i), then one arrives at a finite-difference approximation 
of the eq.(4.2). The proof is finally completed, if one applies classical results on 
numerical solution of the second-kind Fredholm integral equation. 

An example for K (x, t) = A(x — t) , /' (x) = — n ~ / (x) = —nx is 
shown in Fig. 2, where numerical solution is compared with the exact one g (x) = 



8V1 - x 2 (4 + Ax) /(32 + A 2 ) in the case A = 3, (a, 6) = (-1, 1). 
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It should also be noted that some interesting results on the theory and numer- 
ical treatment of hypersingular integral equations have been recently obtained in 
[8-10]. 



5 Some examples from mechanics and acoustics 
5.1 Single crack in linear isotropic elastic medium 

Consider the two-dimensional in-plane problem. Then the displacement field is 
u = (u Xj ,u y ,0), where u x = u x (x,y), u y = u y (x,y) can be represented in the 
Papkovich-Neuber form 



and v — A/2(A + /i) is the Poisson ratio (A and \i are elastic constants). 

Let a constant normal load a yy = (r [r xy = 0) be applied symmetrically to 
internal faces of the straight-line crack y — 0, x G (—a, a). Then, due to a 
symmetry, the problem can be reduced to the upper half-plane y > with the 
boundary conditions on the boundary line y = 

T xy = 0, — oo < x < oo ; u y = 0, \x\ > a ; a yy = —uq, \x\ < a . (5.2) 

Solution of eqs.(5.1) can be constructed with the use of Fourier transform (all 
capital letters are Fourier transforms of corresponding functions): 




where the potentials y? and ip are harmonic 



A<^ = 0, 



= , 



(5.1) 



$ (s, y) = A (s) e-\ s \\ * (s, y) = B (s) , 



where the well known properties 



/(*) 



F (s) , /" (x) 



s 2 F(s), so A(<p,4>) 



dy 2 



($, *) - s 2 ($, 
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have been used. 

The two constants A(s) and B(s) can be determined from the boundary con- 
ditions (5.2) using 

dy 



2-^- + (l-x) — 
dxdy dx 



<Tyy = 2 f jL-^+(l- X )(\ + 2 f J,y 



(5.3) 



It follows from (5.2)i that 



2{-is) {-\s\) A(s) + (1 - x) Hs) B(s) = A(s) 



i-x 

2|s| 



5(a) . (5.4) 



Introducing a new function 



u 



y i 



u y = 



Ixl > a 



,0) = / git) 



e lst dt 



(5.5) 



using eq.(5.4), we obtain the constants A (s) and B (s) in terms of g (x) 



A(s) 



2 (*) «t , 5 (s) 



g (t) e lst dt . 



/From (5.3) one obtains the Fourier transform of the normal stress &yy: 



Zyy (s, 0) = 2pis 2 A(s) + (1 - X )(A + 2 M ) 



1 - i/ 



s| g(t)e tst dt, 



and the boundary condition (5.2) 3 leads thus to the integral equation (|x| < a) 



a . 

J g{t)K{t- x) dt = (1 ~ V) a () . K (.r) = 



s e 



c ds = -— . (5.6) 



x" 



Equation (5.6) is the characteristic hypersingular integral equation 



9 it) 



f g_0 



7T(1 - Z/) 



do, |x| < a 



(5.7) 



(x - ty n 

wit li respect to the function (7 (x) which is, due to eq.(5.5), the relative opening 
of the crack's faces. Solution of eq.(5.7) is (see Section 3) 

g (x) = u y (x, 0) = — (1 — v)\J a? — x 2 , \x\ < a , (5.8) 

which can also be obtained by the proposed direct numerical collocation method 
(see Fig.l). Note that (5.8) coincides with a classical solution of the considered 
problem. 



12 



5.2 Diffraction by a thin rigid screen 

Acoustically hard thin screen is placed horizontally on the interval —a < x < 
a, y = 0. A plane acoustic wave is normally incident on the screen: p tnc = 
e ik y ^ _ ^ w h ere ^ j s the angular wave freguency, c is the wave speed. 
The total pressure is then a sum of the incident and the scattered wave fields 
p (x, y) = p mc Jf-p sc ; where all three components of the wave field satisfy Helmholtz 
equation 

Ap+k 2 p = 0. (5.9) 

By analogy to the previous problem, the Fourier images in the lower (x < 0) 
and upper (x > 0) half-plane are, respectively (7 (s) = \/ s 2 — k 2 ): 

Pi c (s, y) = A (s) ,y<0; P^ c (s,y) = B (s) e"^ , y > . (5.10) 

The coefficients A and B are defined from boundary conditions at y = 

dp + (x, 0) _ dp.jx, 0) _ Q ^ dp™ (x, 0) _ dpt c (x, 0) _ ik ^ ^ 



dy dy dy dy 

dp+ (x, 0) _ dp- (x, 0) dp^ (x, 0) _ dp s J (x, 0) 



\x\ > a (5.12) 



dy dy dy dy 

p+(x,0) =p_(x,0) , p+ (x, 0) = p s l (x, 0) , \x\ > a . (5.13) 

The conditions (5.11), (5.12) imply that dp+/dy = dp s _f/dy for y = and all 
—00 < x < 00, thus 'jA = —'-fB B (s) = —A(s). Now it is obvious from 
(5.10) that pl c (x,0) = — p+(x,0). It thus implies, together with (5.13), that 
pi c ( x , 0) = -p s _£ (x, 0) = 0, \x\ > a. 

Introducing a new function g (x) 



, \x\ > a 
g (x) , \x\ < a , 



p sc (x,0) = 
it follows from (5.10) that 

a 

A(s) = J g(t)e* st dt, 
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and the boundary condition (5.11) leads to the integral equation involving func- 
tion g (x) 

a oo 

J g{t)K(x -t)dt = -ik, \x\ < a , K (x) = — J 7 (s) e~ isx ds , (5.14) 

—a —00 

where the kernel 

00 00 
K( x ) = —J \s\e~ isx ds + — j (Vs 2 - k 2 - \s\) e~ isx ds = 

(5.15) 



1 



-00 
00 



TlX 2 IT 



+ K (x) , K (x) — — / (V s 2 — k 2 — s) cos (sx) ds , 



and the regular kernel K Q (x) is differentiable for all x 7^ 0, with a weak singularity 
in origin: K$ (x) = O (In \x\) , x — * 0, that permits application of the proposed 
numerical method. An example for some values of k is shown in Fig. 3. 



6 Conclusions 

l. The proposed method is an efficient alternative to a standard reduction to 
infinite systems of linear algebraic equations. This is based on direct numerical 
treatment of hypersingular integrals, and reduces the problem to a finite system of 
linear algebraic equations. Its principal merit is that there is no need in numerical 
computations when calculating elements of respective matrix. 

2. In some cases the proposed method permits explicit analytical solution 
of respective hypersingular integral equation. Otherwise, the proposed method 
provides an efficient numerical treatment, and convergence of the algorithm has 
been proved when the step of the mesh tends to zero. 

3. The considered examples show a good precision of the method, since with a 
few decades of nodes on the mesh the obtained numerical results almost coincide 
with respective analytical solution, in the cases when the latter is known. 
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Legends to Figures 

Figure 1. Comparison between exact and numerical solutions of the char- 
acteristic equation with f'(x) = —n ~ f(x) = —nx: - exact solution 
g(x) = y/l — x 2 , numerical solution; n = 40. 



Figure 2. Comparison between exact and numerical solutions of the full equa- 
tion with f'(x) = —n ~ f(x) = — 7T x, K (x,t) = A(x — t), A = 3: 
exact solution, numerical solution; n = 40. 



Figure 3. Comparison between exact and numerical solutions of the equation 
(5.14)-(5.15), (a = 1): 1 - k = 0.5, 2 - k = 1.5, 3 - k = 2.5. The factor 
— nik is omitted. 
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